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We examine spinor Bose-Einstein condensates in optical superlattices theoretically using a Bose- 
Hubbard Hamiltonian that takes spin effects into account. Assuming that a small number of spin- 
1 bosons is loaded in an optical potential, we study single-particle tunneling that occurs when 
one lattice site is ramped up relative to a neighboring site. Spin-dependent effects modify the 
tunneling events in a qualitative and quantitative way. Depending on the asymmetry of the double 
well different types of magnetic order occur, making the system of spin-1 bosons in an optical 
superlattice a model for mesoscopic magnetism. We use a double-well potential as a unit cell for a 
one-dimensional superlattice. Homogeneous and inhomogeneous magnetic fields are applied and the 
effects of the linear and the quadratic Zeeman shifts are examined. We also investigate the bipartite 
entanglement between the sites and construct states of maximal entanglement. The entanglement 
in our system is due to both orbital and spin degrees of freedom. We calculate the contribution of 
orbital and spin entanglement and show that the sum of these two terms gives a lower bound for 
the total entanglement. 



PACS numbers: 03.75.Mn,03.75.Lm,03.75.Gg 



I. INTRODUCTION 

Ultracold atoms can be trapped via the ac Stark effect 
in optical lattices, which are created by counter/propa- 
gating laser beams; in case there are only a few atoms 
per site they build up optical crystals. These systems 
offer a unique combination of experimental and theoret- 
ical accessibilities [l|. They can be manipulated with a 
very high degree of accuracy and versatility so that they 
serve as quantum simulators, i.e. they can be used to 
simulate complex problems in many-body physics. Ul- 
tracold atoms in optical lattices offer robust quantum 
coherence, a unique controllability and powerful read- 
out tools, such as time-of-flight measurements @, H| or 
fluorescence imaging 0]. 

Trapping ultracold atoms in conventional magnetic 
traps leads to frozen spin degrees of freedom such that 
the atoms behave effectively as spinless particles. If the 
atoms are trapped by optical means only, the atoms keep 
the extra spin degree of freedom and the Bose-Einstein 
condensate becomes a spinor condensate. The spinor de- 
gree of freedom on alkaline gases corresponds to the mani- 
fold of degenerate Zeeman hyperfine levels. The ground- 
state properties of spinor Bose-Einstein condensates in 
single traps were investigated in Refs. [!, H|. 

We study the behavior of spin-1 atoms in optical su- 
perlattices, in particular, an optical lattice that is formed 
by overlapping two standing- wave laser fields with a com- 
mensurate wavelength ratio of 2. The resulting lattice is 
an array of optical traps with a double- well structure. We 
model the case when each double-well potential is filled 
with a small number of spin-1 bosons. Spin-1 bosonic 
atoms in a double-well potential can be described by a 
variant of the two-site Bose-Hubbard Hamiltonian Q. 



This model allows examining the interplay between the 
kinetic energy (embodied by the tunneling strength be- 
tween the sites) and the particle interaction (covered by 
the on-site interaction, i.e. the interaction within the 
wells). Furthermore, it is possible to include an energy 
offset between the sites, and the spin-1 Bose-Hubbard 
model additionally contains a term that incorporates 
spin-dependent interactions. This term penalizes high- 
spin configurations on individual lattice sites in the case 
of antiferromagnetic interaction between the atoms (e.g. 
for 23 Na) and low-spin configuration in the case of ferro- 
magnetic interactions (e.g. for 87 Rb). 

The two-site Bose-Hubbard model for spinless bosons 
can be used to describe the transfer of single Cooper pairs 
in small Josephson junctions, i.e. the physics of Cooper- 
pair staircases [a-[l0] . With ultracold atoms in optical su- 
perlattices this model was realized and was shown to give 
rise to a single-atom staircase [TlT - fl5l |. This is achieved 
by monitoring the particle number in either of the wells 
for different values of the energy offset. In the case of 
small tunneling strength, the difference in the number of 
atoms in the two wells does not change smoothly when 
the energy offset is varied but is characterized by a step- 
like behavior. Jumps from one step to the next signal 
the tunneling of a single atom. In this paper, such single- 
atom staircases are examined for spinor condensates. De- 
pending on the energy bias, different types of magnetic 
order occur, and the system of spin-1 bosons in an optical 
superlattice becomes a model for mesoscopic magnetism. 
A specific example of how this mesoscopic magnetism can 
be observed in experiments is presented in Fig. [T] This 
figure shows the difference between bosonic staircases for 
two spin-1 bosons for configurations with different total 
spins. If the total spin is Stot = 2, the spins of the two 
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atoms are parallel and for antiferromagnetic interactions, 
as in the case of 23 Na, being in the same well costs ex- 
tra energy. Therefore the Stot = 2 configuration switches 
later (i.e. at higher energy offset) to the state with both 
atoms in the same well. In the ferromagnetic case (such 
as 87 Rb), the curves for Stot — and Stot = 2 will be 
exchanged. 
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FIG. 1. (Color online) Two spin-1 bosons with antiferromag- 
netic ordering in a double- well potential. Here ur is the occu- 
pation number of the right well, and e characterizes the energy 
offset between the two wells (t/Uo = 0.05 and U2/U0 = 0.1). 
Depending on the total spin of the system, bosonic staircase 
transitions occur at different bias voltages. Note that both 
states with S to t = (red, solid line) and Stot = 2 (blue, 
dashed line) have symmetric orbital wavefunctions with re- 
spect to particle exchange. The difference in the occupation 
numbers arises due to spin-dependent interactions and not 
due to a different orbital symmetry of the states. Thus, a 
measurement of the spin-dependent bosonic staircases pro- 
vides a demonstration of mesoscopic magnetism. 



Spin-1 atoms also allow stronger quantum correlations 
between the wells compared with the case of spinless 
bosons. For spinless bosons it has been noted that parti- 
cle fluctuations between the left and the right well lead to 
entanglement between the wells (see e.g., Refs. 0, [I?} 
and references therein). In addition to this orbital en- 
tanglement, spin-1 atoms allow spinor entanglement. In 
this paper, the quantum correlations between the wells 
are examined for different values of the energy offset and 
different ratios of the tunneling strength relative to the 
on-site interaction. We give a lower bound for the en- 
tanglement between the wells by estimating the amount 
of orbital and spinor entanglement separately. At this 
point, we consider entanglement mainly as a theoretical 
characterization of the many-body state of the system. 

The paper is organized as follows: In Sec. [Til the two- 
site Bose-Hubbard model for spin-1 atoms is introduced 
and is given explicitly for a small number of bosons. In 
Sec. Mil the physics of the bosonic staircases is discussed, 
and in Sec. lIII Cl the effect of magnetic fields is included. 
In Sec. IIV1 the bipartite entanglement for the two-site 
Bose-Hubbard model is examined. The total entangle- 
ment between the sites depends on orbital and spin de- 
grees of freedom. We obtain a lower bound of the total 



entanglement, which is given by the sum of the orbital 
entanglement and the spin entanglement. 



II. TWO-SITE BOSE-HUBBARD 
HAMILTONIAN FOR SPIN-1 ATOMS 

The atoms we have in mind are alkali-metal atoms, 
such as 23 Na and 87 Rb. Degenerate gases of alkali-metal 
atoms are weakly interacting systems, but due to the 
confining lattice of counterpropagating laser beams some 
of the atoms are forced to be very close to each other 
and, thus, to become strongly interacting. Spin-1 bosonic 
atoms in a double-well potential can be described by a 
variant of the two-site Bose-Hubbard Hamiltonian [7lll8j. 

H ° = Y 12 niim-^-tY^iLlRa + RlK) 

i=L,R a 

+e{n L -n R ) + 1 ^ £ (s? - 2n/) , (1) 

i=LM 

where (ij.) and R a (iij.) are annihilation (cre- 
ation) operators for atoms in the hyperfine state a € 
{ — 1,0,1} in the left or right well and i%l = ^ a L\L a 
{ n B = J2a R-o-Rcr) is the atom number at the left (right) 
site. The annihilation and creation operators obey the 
canonical commutation relations [Li, it] = Rlj = 

and [R^L]] = [L h Rt] = 0. § L = L\f aa ,L a , is the 

total spin on the left site and the total spin on the right 
site is Sr — YI0-0-' R0-T0-0-' R<t' , where T aa i are the usual 
spin-1 matrices. 

The on-site repulsive interaction is described by the 
first term in Eq. (JTJ) and parametrized by XJq. The second 
term embodies the spin-symmetric tunneling between the 
wells and t is the hopping matrix element between the 
lattice sites; e characterizes the difference in on-site en- 
ergy between the sites. The term proportional to U2 de- 
scribes spin-dependent interactions: It penalizes nonzero 
spin configurations on individual lattice sites in the case 
of antiferromagnetic interactions (e.g., 23 Na) and favors 
high-spin configurations in the case of ferromagnetic in- 
teractions (e.g., 87 Rb). 

The parameters can be controlled by adjusting the in- 
tensity of the laser beams; it is possible to move from 
regimes of strong tunneling (U$ <C t) to regimes of very 
weak tunneling (t <C Uq). For bulk lattices, it has been 
shown theoretically [l~8| and experimentally that the 
system can be in a Mott-insulating regime (for t <C Uq) 
and in a superfluid phase, where the kinetic energy dom- 
inates (for Uq <C t) , and that it is possible to switch from 
one regime to the other by tuning the laser strength. 

Whereas the ratio of Uo/t can be changed, the ratio 
U0/U2 is fixed for all lattice geometries. £/ 2 is given by the 
difference in the scattering length of two spin-1 bosons in 
the case that their spins couple to the total spin two, and 
the scattering length in the case that their spins couple 
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to the total spin zero. This leads to an estimated ratio 
of[/ 2 /[/ = 0.04for 23 NaQ. 

The Hamiltonian ([TJ conserves the particle number, 
the z projection of the total spin, i.e. (Sl + Sr) z , and 
the total spin (Sl + Sr) 2 . In each well the bosonic angu- 
lar momenta couple, where symmetry constraints require 
rii + Si to be even, and the resulting spins of both wells 
couple to a total angular momentum S tot = Sl + Sr. 



A. Two spin-1 bosons 

Since the hopping term in Eq. ([T]) conserves the abso- 
lute value of the total spin S to t = \Sl+ Sr\, the Hilbert 
space decomposes in orthogonal subspaces that do not 
mix, i.e. 

H = H(S tot = 0) © H(S tot = 1) © H(Stot = 2) . 

In the case of two spin-1 bosons, the Hilbert space is 
seven dimensional 



Ei) 


= |{2,0},{0,0},0), 


E 2 ) 


= |{1,1},{1,1},0), 


E 3 ) 


= |{0,2},{0,0},0), 


Ei) 


= |{1,1},{1,1},1), 


E 5 ) 


= |{2,0},{2,0},2), 


E 6 ) 


= |{1,1},{1,1} ! 2), 


E 7 ) 


= |{0,2},{0,2},2), 



using the notation \{til, ur}, {Sl, Sr}, Stot)- These ba- 
sis vectors belong to three orthogonal subspaces, 



H = {E 1 ,E 2 ,E 3 } 



s tot =o 



{Ek} 

Stat = l 



{E 5 , E e , E 7 } 

S tot =2 



To examine the ground-state properties of this system the 
Hamiltonian needs to be calculated and diagonalized for 
each subspace separately. To calculate the off-diagonal 
elements of the Hamiltonian it is necessary to write the 
elements of the whole system as product of the single- well 
wavefunctions, e.g., 

|{2,0},{2,0},2) = 
\n L = 2, S L = 2 , Si, = 0) ® \n R = 0, Sr = 0, S 2z = 0) 

and 

|{M},{1,1},2) 

= E ^mWm)! 1 ' 1 '™)®! 1 ' 1 '-™) 



m=-l,0,l 
1 



11,1,-1} 

+ 2|1,1,0) <8|1,1,0} + |1,1,-1><8|1,1,1» 



component. The single- well wavefunctions need to be 
written in terms of single-particle creation operators. For 
two spin-1 bosons this can be performed using the stan- 
dard Clebsch-Gordan coefficients. 
The diagonal elements are given by 

(E 1 \H\E 1 )=2e + U -2U 2 , 
(E 3 \H\E 3 ) = -2e + U -2U 2 , 
(E 5 \H\E 5 ) = 2e + U + U 2 , 
(E 7 \H\E 7 )=-2e + U + U 2 , 
(E 2 \H\E 2 ) = (E4H\E 4 ) = (E 6 \H\E 6 ) = . 

Due to the conservation of the total angular momentum 
the Hamiltonian is block diagonal. The only nonvanish- 
ing tunneling elements are 

(Ei\H\E2) = {E 2 \H\E 3 ) = -V2t, 
(E 5 \H\E 6 ) = (E 6 \H\E 7 ) = -yfu. 



B. Higher boson numbers 

The Hilbert space for three spin-1 bosons is given by 
the direct sum of the following subspaces: 



H — {Ei,E 2 ,E 3 , E4, E 5} E 6 } 
®{Eg, E10, En,Ei 2 }. 



{E 7 , E s } 



The subspace {E 1 ,E 2 ,E 3 ,E i ,E 5 ,Eo} belongs to \Sl + 
Sr\ = 1, the subspace {E 7 ,E 8 } belongs to |<Sl + >Sr| = 2, 
and the subspace {Eg, Ei Q , En, E 12 } belongs to \Sl + 
Sr] = 3. The basis vectors are given by 



\Ei) 
\E 2 ) 
\E 3 ) 
\E±) 
\E 5 ) 
\E 6 ) 
\E 7 ) 
\Es) 
\E 9 ) 
\Eio) 
\En) 
\Ei 2 ) 



|{3,0}, 
|{2,1}, 
|{2,1}, 
|{1,2}, 
|{1,2}, 
|{0,3}, 
|{2,1}, 
|{1,2}, 
|{3,0}, 
|{2,1}, 
|{1,2}, 
|{0,3}, 



{1,0},1) 
{2,1},1) 
{0,1},1) 
{1,2},1) 

{1,0}, 1> 
{0,1},1) 
{2,1}, 2) 
{1,2}, 2) 
{3,0}, 3) 
{2,1},3> 
{1,2},3) 
{0,3}, 3) 



Sr}, St, 



where we have chosen the S z = component for con- 
venience because the energy does not depend on the S z 



again using the notation \{ul, tir}, {Sl, 

The Hamiltonian is block diagonal in the basis given 
above as in the case of two bosons. The Hamiltonians be- 
longing to I Sl + Sr I = 2 and | Sl + Sr \ = 3 are quite sim- 
ilar to the spinless case and the case of two spin-1 atoms. 
To calculate the off-diagonal elements of the Hamiltonian 
we need to know how n spin-1 bosons couple to a total 
spin S with a ^-projection S z . For three spin-1 bosons 
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this is performed in the Appendix. The off-diagonal ele- 
ments for \Sl + or | = 2 and \Sl + 6r| = 3 are 

(E 7 \H\E 8 ) = -t, 

(E 9 \H\E W ) = (EulHlEu) = -Vst, 
(E W \H\E U ) = -2t . 

The diagonal elements are given by 

(E 7 \H\E 7 ) = (E W \H\E W ) = U Q + e + U 2 , 
(E S \H\E 8 ) = (E u \H\E n = U -e + U 2 , 
(E 9 \H\E g ) = U + 3e + 3U 2l 
(E 12 \H\E 12 ) = 3U -3e + 3U 2 . 

The Hamiltonian belonging to the Hilbert space \Sl + 
Sr\ = 1 exhibits a richer structure and differs from the 
spinless case. This is because the term —t^^L^R^ + 
W a L a ) describes tunneling between several basis vectors, 
e.g. between \E\) and \E 2 ) as well as between \E\) and 
\Ez). Because the energy does not depend on the S z 
projection we can set S z = 0. The basis vector \E%) is 
given by 

|{3, 0}, {1, 0}, 1> = |3, 1, 0) ® |0, 0, 0> 
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LUL f L\ - 

and the basis vector \E 2 ) is given by 
|{2,1},{2,1},1) 

E E 



10), 



d 1,0) Wl J2,2,m) 

(l,m),(l,n)l ' 1 / 



m=— 2,...,2 n=— 1,0,1 

~3 

To 

~3 

To 



2,1,-1) <g> 1 1, 1, 1> 
2,1,1)® |1, 1,-1) 



12,1,0)0 11,1,0) 



o 



Aft f tot 
jq-^I-^O-"--! 



—V_ X L Q R\ - d —L < _ 1 L[R{ ) 



t ft At 



|0). 



Now, we can calculate the corresponding nondiagonal el- 
ement of the Hamiltonian, 

(£ x |ff|£ 2 ) = -t (Et^&tA* + R\L a )\E 2 ) =~\J\t. 

a 

The basis vector \E%) is given by 
|{2, 1},{0, 1},1)= |2, 0,0)® 1 1,1,0) 



-L^lIrI - -j= (z^ r\ 



10), 



and the corresponding nondiagonal element of the Hamil- 
tonian is 



(£i|#|£i) 



Note that the off-diagonal elements of the Hamiltonian 
depend on the spin configurations, also, the off-diagonal 
elements do not depend on the strength of the spin- 
dependent interactions U 2 . Similar calculations lead to 
the remaining off-diagonal elements, 



(E 2 \H\E 4 ) = -t 

(E 2 \H\E 5 ) = (E 3 \H\E A ) 

2t 
~P 



2y/dt 



(E 3 \H\E 5 



(E4H\E 6 ) = (E 1 \H\E 2 )=- X /-, 



t. 



(E 5 \H\E 6 ) = (E!\H\E 3 ) = 

The diagonal elements are given by 

(E 1 \H\E 1 )=3U +3e-2U 2 , 
(E 2 \H\E 2 ) =U + e + U 2 , 
(E 3 \H\E 3 ) =U + e-2U 2 , 
{E A \H\E A ) = U -e + U 2 , 
(E 5 \H\E 5 ) = U -e-2U 2 , 
(E e \H\E e ) =3U -3e-2U 2 . 

Higher boson numbers lead to analogous expressions 
that are used in the following but are not given here. 



III. BOSONIC STAIRCASES 

The two-site Bose-Hubbard model may be used to 
model Cooper-pair staircases [8HTo| relevant for small 
Josephson junctions. In the case of ultracold atoms the 
same effect gives rise to single-atom staircases 
Here we present such staircases for spin-1 atoms. 



A. General treatment 

The Hilbert spaces decompose into different subspaces 
according to the total spin of the system. Different sub- 
spaces are not mixed by ramping up the energy dif- 
ference between the double wells and behave in a dif- 
ferent way. In the case of two bosons this is shown 
in Fig. [T] The different widths of the steps centered 
at e = correspond to different values of total spins 
per site. At e = 0.5, state |{1, 1}, {1, 1}, 2) is ener- 
getically lower than state |{0, 2}, {0, 2}, 2), which makes 
the step broader for Stot = 2. On the contrary, at 
e = 0.5, state |{1, 1}, {1, 1}, 0) is energetically higher than 
state |{0, 2}, {0, 0}, 0), which makes the step narrower for 
S to t = 0. 

In general, depending on the sign of U 2 , states with 
high single-well angular momenta get penalized or get 
favored. If U 2 > (such as, e.g., for 23 Na), nonzero 



spin configurations get penalized. In the case of Rb, 
U2 is negative and spin-dependent interactions lead to 
the opposite effect: High-spin configurations are favored, 
and the corresponding steps are broader. Therefore, in 
the ferromagnetic case, the curves for Stot = and Stot 
in Fig. [T] will be exchanged. 

The exact position of the steps can be calculated in 
the atomic limit, i.e., t = 0. The step positions generally 
depend linearly on t/2. For some spin configurations, 
e.g., odd atom number, lowest possible total spin, and 
antiferromagnetic interactions, the step positions do not 
depend on spin-dependent interactions. 

For higher boson numbers, the richer structure of the 
off-diagonal elements means that the variance in the par- 
ticle number depends on the total spin and the energy 
offset. In the case of three bosons (Fig. [2j, the step at 
e = is not shifted due to symmetry reasons, whereas 
the steps at e = 1 and e = — 1 are shifted linearly. At 
the same time, the steps belonging to Stot = 3 are not as 
sharp as the steps belonging to S to t = 1, i-e., the curve 
of the variance in is broader in the case of S to t — 3. 

<« 2 > 




FIG. 2. (Color online) Bosonic staircase for three spin-1 
bosons with antiferromagnetic ordering in a double-well po- 
tential (t/Uo = 0.05 and U2/U0 = 0.1): Stot = 1 (red, solid 
line) and Stot = 3 (blue, dashed line). (Inset) Variance in the 
particle number in the left well for the step around e = 0. 



The staircases for different total spins may be used to 
arrange spin-1 atoms in a two-dimensional (2D) super- 
lattice according to their spin degrees of freedom (see 
Fig. El). 
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FIG. 3. (Color online) A possible way to separate the S to t = 2 
spin component from the Stot = spin component in the case 
of antiferromagnetic interactions between the atoms: The po- 
tentials in the x and y directions are manipulated separately. 
In the first step (b), the energy offset between the double 
wells is lifted until the bosons combining with the total spin 
Stot = separate while the bosons belonging to Stot = 2 still 
remain in the same site, (c) Next, the wells are separated 
by a large potential barrier and tunneling is suppressed, (d) 
An additional laser is switched on and the bosons coupling to 
Stot = 2 distribute in the resulting double well. The switch- 
ing is assumed to happen adiabatically such that the system 
can be regarded to be in the ground state at every instant. 




-1.0 -0.5 0.0 0.5 1.0 



FIG. 4. (Color online) Energy spectrum of two spin-1 bosons 
in a double well with strong tunneling (t/Uo = 0.5 and 
U2/U0 = 0.1): Stot = subspace (red, solid lines), Stot = 2 
(blue, dashed lines), and Stot = 1 (dotted line). 



B. Beyond ground-state analysis 

The gap between the ground state and the first excited 
state in the energy spectrum depends strongly on the 
tunneling between the sites (see Figs. 2] and EJ. For finite 
temperatures the density matrix describing the system, 
thus, is highly mixed for small tunneling parameters, and 
the ground-state behavior only dominates if tunneling is 
sufficiently strong. 



C. Magnetic field included 

The effect of a magnetic field can be included in the 
model |T]) by adding a term to the Hamiltonian thath de- 
scribes the coupling of the spins to the magnetic field ■ 
The first contribution of a magnetic field B = (0, 0, B) is 
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-1.0 -0.5 0.0 0.5 1.0 



e/U 

FIG. 5. (Color online) Energy spectrum of two spin-1 bosons 
in a double well with weak tunneling (t/Uo = 0.05 and 
U2/U0 = 0.1). Color code as in Fig.g] 



P/U 




FIG. 7. (Color online) Critical magnetic field p — q\ibB 
above which the staircase for two bosons shows a discontin- 
uous behavior signifying spin- flip transitions (U2/U0 = 0.1). 



a regular Zeeman shift in the energy levels: 

H = H Q +p m ™ fl ™ = H n+P S T 

i=L,R <J 

where p = gusB and h ia is the particle number opera- 
tor for the ith site that gives the number of bosons in the 
mth hyperfine state. The linear Zeeman shift changes 
the overall state considerably. The energy eigenvalues 
belonging to Stot 7^ split into multiplets because the 
hyperfine levels are no longer degenerate (see Fig. [5]). For 



: ^sj!^ 









-1.0 -0.5 0.0 0.5 l.i 



e/U 

FIG. 6. (Color online) Linear Zeeman shift of the energy levels 
in the energy spectrum of two spin-1 bosons (B/Uo = 0.05, 
t/Uo = 0.05 and U 2 /U = 0.1). The energy levels of Fig. g] 
split into spin multiplets. The red arrows denote ground state 
level crossings. Color code as in Fig. [U 

a given tunneling strength, there is a critical magnetic- 
field strength that leads to ground-state level crossings. 
Such level crossings correspond to spin-flip transitions, 
i.e., the ground-state energy is continuous, but the ex- 
pectation values of the particle number of single sites 
and of the magnetization are not. This means that the 
overall ground state of the system does not belong to the 
same z projection of the total spin for all values of the 
energy offset e. Figure [7] shows the critical value of the 
magnetic field in the case of two bosons. 



However, spin-non-conserving collisions are negligible 
over the lifetime of the condensate, and the total mag- 
netization is a conserved quantity on the time scale of 
the experiment [111, [22[ • For a given magnetization the 
properties of the system are not altered by the linear Zee- 
man effect; the whole spectrum merely is shifted. Only if 
one is interested in comparing different magnetizations, 
the linear Zeeman effect has to be taken into account. 
In a series of experiments with a given magnetization, it 
is necessary, therefore, to include higher-order contribu- 
tions in the magnetic field. The quadratic Zeeman effect 
arises because the hyperfine spins characterizing ultra- 
cold atoms are mixtures of electron and nuclear spins. 
Since the magnetic field approximately couples only to 
the electron spin, the Zeeman effect is nonlinear in the 
field but , typically, can be described by a sum of linear 
and quadratic terms. 

For each of the subspaces belonging to different mag- 
netizations S 1 ' ', there is a separate effective Hamiltonian 

H q =H Q + q TO -"- • ( 2 ) 

i=L,R <t 

The magnitude of the quadratic Zeeman shift is given by 
q = q B 2 , where e.g. q = h x 390 Hz/G 2 for Na HJ. 

In the case of two bosons the system with magnetiza- 
tion S z ot — possesses the most interesting structure 
because the Hilbert space is composed of states with 
different total spins. For S z ot — 2, the quadratic Zee- 
man shift does not alter the staircase since it leads to 
a homogeneous shift in all the energy levels. The stair- 
cases at different magnetic fields are shown in Fig.[5J For 
S ° = 0, the step positions depend in a nonlinear way 
on the magnetic field strength. It is no longer possible 
to read them off in the atomic limit (i.e., t = 0) because 
the existence of the quadratic Zeeman shift leads to ad- 
ditional nondiagonal elements in the Hamiltonian. Note 
that the quadratic Zeeman effect does not eliminate the 
difference in the two staircases, which is the main mani- 
festation of mesoscopic magnetism. 
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FIG. 8. (Color online) Two spin-1 bosons with antiferromag- 
netic ordering in a double-well potential (t/Uo = 0.05 and 
U 2 /U = 0.1). (Dashed lines) 5"*°' = for different mag- 
netic fields q = qoB 2 (short dashes q/£/o=0.2, long dashes 
q/f7o=0.f and medium sized dashes q/Uo=0). (Solid line) 
gtot _ 2 j n ^ n j s case tn e staircase does not depend on the 
magnetic field. The difference in this staircase to the ones 
with S\ ot — 0, which is the main manifestation of mesoscopic 
magnetism, persists in the presence of the quadratic Zeeman 
effect. 



Due to the fact that the quadratic Zeeman shift does 
not commute with the operator of the total spin S tot , 
the eigenstates of the Hamiltonian given in Eq. (|2) arc 
no longer eigenstates of S to t- For B ^ 0, the ground state 
of the system is a superposition of different eigenstates of 
Stot, i-e., states with different Stot hybridize (see Fig. [9]). 
For certain values of the energy offset e (e.g., c/Uq = 1 
and e/Uo = — 1 for four bosons), the appearance of a 
magnetic field changes the ground state strongly. This 
reflects the specific spin configurations. 




FIG. 9. (Color online) Expectation value {Stot) of the system 
for four bosons for different magnetic fields q — qoB 2 (5' oi = 
0, t/U = 0.05 and U 2 /U = 0.1). 

The quadratic Zeeman shift also changes the overall 
spectrum for a given magnetization qualitatively such 



that, in the case of thermal occupation of the double 
well the density matrix of the system changes consider- 
ably. For q — 0, the ground state nearly is degenerate 
with the first excited state, whereas the gap widens for 
finite values of q. 

Additionally one can include inhomogeneous magnetic 
fields. 

H = H q + AB (S Lz - S Rz ) , 

where AB describes the strength of the field gradient. 
The magnetic-field offset AB changes the Hamiltonian if 
Stot z 7^ 0. For some configurations, e.g. two bosons in 
a double well, AB merely leads to an overall shift in e, 
i.e., an inhomogeneous magnetic field is equivalent to an 
energy offset e. In general, this is not the case, and AB 
is an additional tool to reshape the staircases depending 
on the spin configuration of the system. 

IV. ENTANGLEMENT FOR SPIN-1 BOSONS 

Entanglement is a unique feature of quantum- 
mechanical systems. Understanding entanglement deep- 
ens our understanding of quantum mechanics and, there- 
fore, is of fundamental interest. Moreover, entanglement 
is a resource for quantum computation and correlates 
separated systems stronger than all classical correlations 
can do. For bipartite pure states entanglement is well 
understood and the different entanglement measures are 
equivalent. In the following we will use the entanglement 
of formation (EOF) [23| as an entanglement measure. 
The EOF is the number of Einstein- Podolsky- Rosen pairs 
asymptotically required to prepare a given state by local 
operations and classical communication. The entangle- 
ment of formation between two qudits (D-dimensional 
objects) in a pure state is given by the von Neumann 
entr opy of the reduced density matrix of each single qu- 
dit [2J]. To calculate the EOF, consider a system con- 
sisting of two parts labeled A and B. Any pure state |\&) 
of the system can be written in the Schmidt decomposi- 
tion [25[ 

D 

|*)=]T Q |^>^f}, 

i 

where {^i 4 , • • • , V'dI an d {V'f 1 ■ ■ - i V'Bl are complete sets 
of orthonormal states of the respective subsystems. The 
system is entangled iff there is more than one nonvan- 
ishing coefficient These coefficients Ci are positive, 
unique and invariant under local operations and, there- 
fore, can be used to quantify the entanglement between A 
and B. The von Neumann entropy of the reduced density 
matrix of each single qudit is given by 

= 5(Tr B |*)(*|) = SC& A |*)<*|) 

D 
1 
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where S indicates the entropy. It ranges from zero to 
log 2 D. The entanglement of formation of two qudits 
with D > 2 thereby exceeds the entanglement of forma- 
tion of two qubits, i.e. higher-dimensional objects con- 
tain more entanglement and violate all Clauser-Hornc- 
Shimony-Holt-inequalities more strongly than qubits. 

In the following we calculate the EOF in a double well 
and examine how much the two sites are entangled. At 
this point, we consider the EOF mainly as a theoretical 
characterization of the many-body state of the system. 

The EOF for typical parameters is shown in Fig. [TO] 
The maximal entanglement exceeds the maximal entan- 



EOF 
2.5 r 





FIG. 11. (Color online) EOF between two wells for four parti- 
cles in the presence of a magnetic field (Sl ot = 0, t/Uo = 0.05 
and U2/U0 = 0.1). For q > 0, even small fields will elimi- 
nate the contribution of the spin degrees of freedom to the 
entanglement. 



FIG. 10. (Color online) EOF between two wells for two bosons 
with antiferromagnetic interactions (t/Uo = 0.1 and U2/U0 = 
0.1) for the total spin S to t = 0. 



glement between two qutrits of log 2 3 ~ 1.585. This is 
due to particle fluctuations. The total amount of entan- 
glement stems from orbital and spin degrees of freedom. 

Magnetic fields have a strong effect on the entangle- 
ment of formation. Figure [11] shows the EOF of four 
bosons at S z ot — 0. For q > 0, the contribution of the 
spin degrees of freedom to the entanglement of formation 
is suppressed already by small magnetic fields. This is 
somewhat surprising, because the system is constrained 
to S z ot , i.e. the state with the strongest spin entangle- 
ment of all states with a given total spin. For q < 0, 
this contribution initially is reduced but then remains 
constant as a function of q. 



A. Two spin-1 bosons 

For two bosons and in the case of Stot — a possible 
orthonormal basis is given by 

{^1,^2,-03} 

= {|{2,0},{0,0},0),|{1,1},{1,1},0),|{0,2},{0,0},0)} 

using the notation \{til, tir}, {Sl, Sr}, Stot}- The de- 
composition, 



is not a Schmidt decomposition, because the vector ^2) 
is a superposition of orthonormal states, 

|{1,1},{1,1},0) = --5=|1, 1,0; 1,1,0) 

1,1;1,1,-1> + -J=|1, 1,-1; 1,1,1), 

using the notation Sl, Sl z ; ur, Sr, Sr z ). The entan- 
glement of formation of \tp2) is given by 



The EOF of |*) is given by 
E(\*)) 

= ~ c i !og 2 c\ - c\ log 2 C3 - 3 ( c 2 ^= I log 2 ( c 2 



V3 



'V3 



= -Y. C * l0 §2 C ? + C 2 l0 §2 3 

i 

3 3 

i i 
-^orbital H~ -^spin- 



(3) 



The total entanglement between the left and the right 
well decomposes in an orbital part and a spin part. The 
orbital part stems from the coefficients that distinguish 
different orbital wave functions. The spin part origi- 
nates from the EOF of the individual basis vectors, each 
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weighted with the coefficient cf. The coefficients Cj de- 
pend on the tunneling strength t, on the on-site inter- 
action Uo, on the spin-dependent interaction U%, and on 
the energy offset e. 

EOF 



2.5 
2.0 




1.5 




— h& 
0.5 





EOF 




FIG. 13. (Color online) EOF between two symmetric wells 
(e = 0) for different values of 11%, i.e., different spin interac- 
tions (solid line U2/U0 = 0.2, dotted line U2/U0 = 0.01, and 
dashed line U 2 /U Q = -0.2). 



FIG. 12. (Color online) EOF between two wells for two bosons 
with very weak tunneling and antiferromagnetic interactions 
(t/U = 0.001 and U 2 /U = 0.1) for the total spin Stot = 0. 

In the limit of weak tunneling t <C Uo the Hamiltonian 
is diagonal in the basis {ipi, ip2-, ^3} an d the ground state 
of a symmetric double-well potential (i.e. e — 0) is given 
by 



EOF 



|*o) = |V'2> = -^=|l,l,0; 1,1,0) 
+ ^|1,1,1;1,1,-1) + -^=|1,1, 



-i; M,i) 



(4) 



which leads to an entanglement of E{\4>2)) — log 2 3, see 

Fig.ru 

In the limit of strong tunneling (i.e. Uq -C t), the 
ground state of the system is 

|* ) = i|^) + -L|^) + i|^). (5) 

For this state the orbital entanglement is given by 
Krbitai = -2j log 2 I - 5 log 2 | = 3/2 and the spin en- 
tanglement is given by E sp i n = (yjl l°g2 3 ~ 0.792. 

Therefore, the total entanglement is B(|* 2 )) ~ 2.292. 
This is not the maximum amount of entanglement that 
can be obtained for this system (see Fig. [T3|) . The max- 
imal entanglement is not the sum of the maximal qutrit 
entanglement and the maximal orbital entanglement be- 
cause the orbital motion leads to particle number fluctu- 
ations and reduces the spin entanglement (see Fig. [Til) . 
The maximal orbital entanglement is realized in the limit 
of strong tunneling; the maximal spin entanglement cor- 
responds to the maximally localized state, i.e., 1^2) • 



B. Three spin-1 bosons 




Em 



10 



t/U 



FIG. 14. (Color online) EOF, -Eapin and -EorMtai between two 
symmetric wells for U2/U0 = 0.1. 



contrast to the case of two bosons, in the weak-tunneling 
case the system is not entangled for large intervals of the 
energy offset e. 




e/U 



FIG. 15. (Color online) EOF between two wells for three 
bosons with antiferromagnetic interactions (t/Uo = 0.1 and 
U2/U0 = 0.1) for the total spin S to t = 1. 



The EOF between the two sites is presented in Fig. [T5] 
for t/Uo = 0.1 and in Fig. [T5]for very weak tunneling. In 



To quantify this effect, we analyze the EOF again in 
detail. The spins Sl and Sr couple to a total spin, for 
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FIG. 16. (Color online) EOF between two wells for three 
bosons with antiferromagnetic interactions (t/Uo = 0.005 and 
U2/U0 = 0.1) for the total spin Stot = !• 



which three absolute values are possible, Stot & {1,2,3}. 

It is obvious that Eq. ^ is applicable for Stot — 2 and 
Stot — 3. The interesting case is Stot = 1< Each state 
with quantum number Stot = 1 can be written as 

l*) = 5>hfc). (6) 

i 

Only two of the basis vectors contain true spin entangle- 
ment, 

|^ 2 ) = |((2,l),(2,l),l))=a|2,2,0; 1,1,0) 

+ /3|2,2,1;1,1,-1)+ 7 |2,2,-1;1,1,1) 

|^> = |((l,2),(l,2),l)) = a|l,l,0;2,2,0> 

+ /3|1, 1,1; 2, 2,-1) +711,1,-152, 2,1) 

with a = — y| and fi = 7 = y^jp Any superposition of 
^2) an d ^3) can be written as 

c 2 |-02>+c 3 |v>3) 

= y/4o? + 4\L)® 11,1,0)* 

+ c 2 /3|2, 2, l;l,l,-l) + c 2 7|2 ) 2,-1; 1,1,1) (7) 

where |L) is the normalized function 
l/Vcla 2 + eg (c 2 a|2,2,0) L + c 3 |2,0,0) L ), which is 
orthogonal to the other vectors appearing in Eqs. 
and (j6|). Therefore, the decomposition Eq. is a 
Schmidt decomposition and the full entanglement of 
formation of \ty) can be calculated, 

E(\V)) = -c? log 2 c? - (c 2 a 2 + c 2 ) l 0g2 (c 2 a 2 + c 2 ) 
-(c 2 /3 2 )log 2 ( C ^ 2 )-(c 2 7 2 )log 2 (c 2 7 2 ) 

- (c 2 a 2 + c 2 ) log 2 (c 2 a 2 + c 2 ) - (c 2 /3 2 ) log 2 (c 2 /? 2 ) 

- (c 2 7 2 )log 2 (c 2 7 2 )-c 2 log 2C 2 . (8) 

It is possible to decompose the entanglement into differ- 
ent contributions and to generalize the expressions for 



-^orbital and -Espin in Eq. © . To calculate the orbital en- 
tanglement we construct the orbital wave function and 
use this to get the EOF of the reduced density matrix. 
The orbital wave function is 

I ^/orbital = 

Cl |3,0) + ^4 + 4\2,l) + v /c 2 + c 2 |l,2)+c 6 |0,3), 

where the quantum numbers refer to {til^r}- So the 
orbital entanglement of formation between the left and 
the right well is given by 

^orbital =-c\ l0g 2 c\ - (c2j + C§) log 2 (cj + cj) 

~ (c 2 + c 2 ) log 2 (c 2 + cl) - cl log 2 cl (9) 
The spin wave function is given by 

Wspin = yfi% + 4\{l, 0}, 1) + c 2 |{2, 1}, 1) 

+ V / ^ _ * r cil{0, 1}, 1) + C4 |{1, 2>, 1), (10) 

where the quantum numbers refer to \{Sl, Sr}, Stot)- 
The EOF of these orthonormal basis vectors is 
£(|{1,0},1)) = £(|{0,1},1)) = and £(|{2, 1}, 1)) = 
£(|{1,2},1)) = -a 2 log 2 a 2 - /? 2 log 2 /? 2 - 7 2 log 27 2 . 
Therefore, the EOF due to spin entanglement is 

E spin = c 2 2 E(\^ 2 )) + clE(\i; 4 )). (11) 

Note that 

(c\o? + cf ) log 2 (c 2 a 2 + 4) + (c 2 2 (3 2 ) log 2 (c 2 /? 2 ) 
+ {cW) log 2 {cW) 

< cl log 2 (c| + cl) + (c 2 ) log 2 ^1 + ^ 

+ c\ log 2 c\~c\E{\^ 2 )) 

= (cl + 4) log 2 (c\ + 4) - 4EQih)), (12) 

where log(l + dz) < dlog(l + z) for d > 1 and z > 
has been used. Because of Eq. (ff"2"|) the entanglement of 
formation is bounded from below (see Fig. [T7)) . 

> ^orbital + E spin (13) 



C. Arbitrary number of bosons 

Let f be a wave function that describes the state of 
N bosons. This wave function can be written in terms 
of a basis, which is ordered according to the occupation 
numbers Nl and Nr, the spin in the left well Sl and the 
spin in the right well Sr, and the total spin Stot', 

D 

i*>-E c «i^)' ( 14 ) 

n=l 
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FIG. 17. (Color online) EOF, E spin , £ orbi tai and E spin + 
Sorbitai between two symmetric wells for three bosons with 
antiferromagnetic interactions (U2/U0 = 0.1) in a symmetric 
double- well potential (e/Uo = 0). 



where J2 n c n = 1 an< ^ D > N is the dimension of the 
basis. We can rearrange this sum by sorting it according 
to the occupation numbers, 

JV 

i*> = E i^>. 

m— 

where \ip m ) is the part of the wave function belonging 
to Nl — m. If N(m) is the number of basis vectors 
belonging to Nl = m, \ip m ) is given by 

N(m) 

\i>m) = E C ( rn )i\ ( t > ( m )i)' ( 15 ) 
i=l 

where c ( m )l — 1 an d c{m)i denote the coeffi- 

cients Ci that belong to Nl = m. Now it is possible 
to generalize Eqs. © and (fTTj) and to define the orbital 
EOF, 

AT N(m) N(m) 
^orbital = -^E °( m ^ lo & E C(m)? (16) 



m— a—l 



and the spin EOF 



(17) 



It is not necessary to to specify which basis vectors 
in Eq. (|14[) belong to which angular momentum con- 
figuration, such as in Eq. ([TUf because the total spin- 
entanglement entropy can be written as a sum over all 
basis vectors. 

In this section we prove that Eq. (|13[) is true for any 
number of bosons in a double well: 



> -^orbital + -^spin 



(18) 



E(\^)) decomposes in a sum over m: -Ed*)) = 



Em-^d^'"})- It is possible to write down the EOF for 



each \ip m ) in the following way: 

E(\Tpm)) = (19) 

which defines a basis for each vector c/)(m)i, 

\<P{ m )i) = E E a (™)ifc a (™)«|ife> O 
ft 1 

= ^a{m) lJ \L,R) j , 

where ■ afj = 1 . To prove Eq. (|18[) for any number of 
bosons, it is necessary and sufficient to show that Eq. (p2 
is true for each E(\ip m )), i.e., 



e(e^) 1o & (ew) 

^ (E c ?) ^ (e c * ) + E c * E 4 «? 

The term on the left-hand side can be rearranged, 

E (e^)i°s 2 (e^) 

= EE(4c?)i°g 2 (4 c <) 



(20) 



j » 



EE K c Di°g2 



T a" c 



2 . c 2 

"J 



2 2 



(21) 



as well as the term on the right-hand side, 

(E lo 8'2 (E + E c ' E a « lo §2 a 

= (E c i) lo g2 fe c ')-E c i lo g2^ 



+ EE(4 c ?) lo g2 K-c 2 ' 

i 3 

Note that, due to Jensen's inequality. 
E ( a ij°i) lo & 



a" c 



2 c 2 

"J 



21 



< c 2 log 2 



: log 2 (E^) 



q iog 2 cf, 

Eq. is fulfilled and, therefore, Eq. CO 
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D. Comparison with the entanglement of particles 

The amount of entanglement shared between two par- 
ties might be lowered by superselection rules [26[ . In case 
two parties share N particles and a particle superselec- 
tion rule applies, the extractable bipartite entanglement, 
i.e., the degree of entanglement one can entangle two ini- 
tially not entangled quantum registers located at A and 
B, is given by the entanglement of particles [jj , 



E P (\y AB )) = Y,Pn E(\* 



ABi 



where \^2bI * s \^Ab) projected onto the subspace of 
fixed local particle number, i.e., n particles for one party 
and n — 1 for the other. 

The entanglement of particles for two bosons in a dou- 
ble well is given by -E S pin in Eq. ([3]) . For three bosons the 
case S to t = 2 and S to t = 3 is trivial, but the case S to t = 1 
is more interesting and is examined. To calculate Ep, we 
write down the projection obeying local particle supers- 
election rules. The projections onto nj, = 3 and m = 
are trivial and do not contribute to Ep. The projection 
onto n L = 2 leads to Eq. © with P 2 = c\ + c\. The 
entanglement contained in this state is given by 



E 



(2)« 
LB i 



C 2 +C 2 



log. 



cW + ci 



¥ 2 



Co + c: 



log 2 



IP 2 



r 2 -v 2 



C 9 + C 3 



loe 



c 27 



and thereby contributes 

P 2 E\^ 2) R ) = - (cja 2 + 4) log 2 (c|a 2 + cf) 

- {clP 2 ) log 2 (cj(3 2 ) - (c 2 2l 2 ) log 2 {cW) 
+ (4 + 4) 1°§2 (4 + 4) 

shows that the equa- 



to Ep. A comparison with Eq 
tion, 



Em) = e 



orbital 



E f 



(22) 



holds for three bosons. This equation also is true for 
higher boson numbers. The contribution of the state (fT5)) 
to Ep is given by 



P m EVS 



(m)v 
LR / 



]c(m)i log 2 



c(m) 



A comparison with Eq. (fTB| shows that Eq. (|2^|) indeed 
holds for all boson numbers. 

The necessity to take a superselection rule into account 
may arise due to several reasons. In some cases the phase 
between states with different local particle occupation 



numbers is not well defined [28]. Consider the bipartite 
state 



\i>e 



AB 



(|l,0) + e*|0,l)) 



(23) 



In case there is no shared reference frame and no tunnel- 
ing between the two parties the phase is not accessible 
experimentally and the state is indistinguishable from an 
incoherent mixture, 



PAB = I (\1,0) (I, 0\ 



|o,i)(o,i|) 



(24) 



Whenever one is concerned with the occupation number 
of massive particles, the detailed properties of the system 
determine which local operations and classical communi- 
cation (LOCCs) are allowed: If tunneling is forbidden, 
LOCCs will conserve the local particle number. In this 
case a local particle number superselection rule must be 
taken into account. A more trivial example is the case of 
a superselection rule for the total particle number |29| . 

In our model ((T|) the phase is well defined due to the 
finite-tunneling amplitude. The amount of orbital entan- 
glement -Eorbitai depends directly on the particle fluctua- 
tions caused by the tunneling between the sites. In the 
absence of tunneling, the orbital entanglement vanishes 
and the superselection rule for the local particle number 
is effectively enforced. 



E. Creation of entanglement structures 

In the case of two spin-1 bosons in a double well the 
state of total spin zero (S to t = \Sl + Sr\ = 0) is sin- 
gled out. First, it can be separated from the Stot — 2 
state due to a different particle distribution within the 
double well in the vicinity of the single-particle tunnel- 
ing resonance (i.e. e/Uo — 0.5). Second, it represents the 
two-qutrit singlet state and thereby contains the maxi- 
mal qutrit entanglement of log 2 3. This distinguishes the 
qutrit entanglement from the qubit entanglement where 
the singlet state and the triplet (Stot)z = state contain 
the same amount of entanglement. 

This can be used to create specific entanglement struc- 
tures in 2D optical superlattices (see Fig. [3]). 



V. CONCLUSION 

We have analyzed the two-site Bose-Hubbard model 
for spin-1 atoms explicitly for small numbers of bosons. 
Starting from the explicit form of the Hamiltonian, we 
have discussed the physics of the bosonic staircases. We 
also have studied the effect of magnetic fields. In the 
following we have examined the bipartite entanglement 
for the two-site Bose-Hubbard model. We have analyzed 
the contribution of orbital and spin degrees of freedom 
and have derived a lower bound of the total entangle- 
ment, which is the sum of the orbital entanglement and 
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the spin entanglement. We compared the entanglement 
of particles and thereby elucidated the meaning of orbital 
entanglement and of superselection rules for the local par- 
ticle number. 

The staircases for different total spins establish a cor- 
respondence of the spatial motion and the spin configu- 
ration. Because the Stot — singlet state of two bosons 
does contain more entanglement than the other eigen- 
states of the system this correspondence can be used to 
construct an entanglement witness in the system: In case 
one detects the typical spatial behavior of the S tot = 
state, one can conclude to have its entanglement. With 
the help of fluorescence imaging, it is also possible to de- 
populate doubly occupied sites in the lattice and thereby 
to build a spin filter. 

We have discussed entanglement between the sites, not 
the entanglement between the individual atoms. Even 
for an occupancy of one, i.e. one atom per site, these are 
different quantities, because the bosons are indistinguish- 
able. Recently it was proposed to measure the entangle- 
ment between (spinless) bosons in an optical lattice [29| 
by standard time-of-flight measurements. Such measure- 
ments do not preserve the information about the entan- 
glement between the individual sites. There are other 
possibilities for examining these systems experimentally. 
First, it is possible to estimate the entanglement by mea- 
surements of the atom positions because these correspond 
to specific spin configurations as we have demonstrated. 
These atom positions can be determined by standard 
time-of-flight measurements or direct fluorescence detec- 
tion of individual sites [||. Furthermore, it is possible to 
detect the spin configurations directly in a nondemolish- 
ing way with the help of the quantum Faraday effect [30l ] . 
Furthermore, it may be possible to relate the entangle- 
ment to additional observable experimental quantities, 
such as magnetization fluctuations in one of the wells, 
in analog y to what has been discussed for noninteracting 
particles [3l| . We plan to explore this question in future 
papers. 
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Appendix A: Coupling of n spin-1 bosons 

The coupling of two spin-1 atoms is given in standard 
textbooks. The coupling of n spin-1 atoms to a total spin 
S with a z-projection S z in terms of n single spins Si Z is 
calculated by the diagonalization of S 2 . 

For three spin-1 atoms we give the connection of the 
basis vectors ordered according to Si Z and the basis vec- 
tors ordered according to S and S z , 
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= 3, S z 


= 3) 
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